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ABSTRACT 

Radial velocity identification of extrasolar planets has historically been dom- 
inated by optical surveys. Interest in expanding exoplanet searches to M dwarfs 
and young stars, however, has motivated a push to improve the precision of near 
infrared radial velocity techniques. We present our methodology for achieving 
58 m s _1 precision in the K band on the MO dwarf GJ 281 using the CSHELL 
spectrograph at the 3-meter NASA IRTF. We also demonstrate our ability to 
recover the known 4 M JUP exoplanet Gl 86 b and discuss the implications for 
success in detecting planets around 1 — 3 Myr old T Tauri stars. 

Subject headings: techniques: radial velocities — planets and satellites: detection 
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Introduction 



Since 1989, astronomers have discovered over 500 extrasolar planets, ~92% of which 



were identified using radial velocity (RV) t echnique; 



improved velocity precision (~1 m s 1 e.g., 



Mayor et al. 



ncreased time baselines and 



20031 ) have led to an abundance 



of ground-based RV planet disc overies including th e lightest known planet around a main 



sequence star (msini = 1.7 M a 



Mayor et al. 



200 9.1) and th e recent claim of a ~3 M e planet 



orbiting within its star's habitable zone (jVogt et al.l 120101 ) . Ongoing exoplanet surveys 
typically focus on slow-rotating FGK dwarfs. The numerous, narrow atomic lines in the 
photospheres of these stars allow for precise RV determination. The SEDs of these stars also 
peak in the optical region of the spectrum, facilitating their study with optical detectors, a 
more mature technology than is available in other wavelength regions, such as near-infrared 
(NIR). 

The past few years, however, have seen a surge of interest in applying precision RV 
techniques in the NIR, in part to explore M dwarfs as potential planet hosts. These late-type 
stars are faint at optical wavelengths but comparatively bright at 1 — 2 /im They are also 
the most numerous stars in our Galaxy; by focusing on FGK dwarfs in most surveys to date 
the majority of stars have been neglected. Also, for a given orbital period and planet mass, 

2/3 

the RV amplitude scales with host mass as M* . Therefore, planets should be easier 
to detect around lower mass M and L dwarfs. Furthermore, including the low mass stars 



expands the available parameter space to include lowe r mass planets. T 



re M dwarfs are 



also interesting targets to the astrobiology community (iTarter et al.l 120071 ). Because of their 
lower luminosity, the habitable zones (HZ) are significantly closer to the central star than 
for higher mass stars and the RV amplitudes of HZ planets are correspondingl y greater, 



thereby increasing the possibility of discovering a habitable Earth-mass planet ( TVogt et al. 



x Data from The Extrasolar Planet Encyclopedia (www.exoplanet.eu) 
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20101). 



M dwarfs may also serve as testing grounds for models of planet formation. In spite 
of all the planet discoveries in the past twenty y ears, a unified mode l of planet formation 



is still elusive. In the core accretion model (e.g., 



Pollack et al. 



19961 ). a planetary core is 



built up through t he accretion of i c es, du st particles, and eventually planetesimals in the 



circumstellar disk. 



Laughlin et al 



( 120041 ) present calculations which indicate that giant 
planet formation around M dwarfs under a core-accretion paradigm is inhibited as a result 
of lower surface densities in the disk beyond the snowline, where Jupiter-mass worlds 
are most likely to form. Thes e planets wou l d thu s take much longer to accrete than the 



typical disk lifetime; however, 



Kornet et al. 



( 120061 ) challenge this result and argue that the 



19971 ) whe rein instabi 



Boss! (120061 ) 



ities 



probability of planet formation actually increases with decreasing stellar mass as a result of 
a more efficient particle redistribution around lower mass s tars. Anot her avenue for giant 
planet formation is the gravitational instability model (e.g.,|B_< 
in the circumstellar disk can lead to fragmentation and eventual collapse, 
argues that gravitational instabilities can produce Jovian worlds around low-mass stars very 
rapidly. Current optical RV surveys of M dwarfs support t he hypothesis th a t low-mass star s 



are le ss likely to harbor planets than Sun-like stars (e.g., 



Endl et al. 



2006 



Johnson et al. 



20101 ) . Johnson et al. report that the occurrence rate of giant planets around M dwarfs 
msini > 0.3 M JUP , ji < 2.5 AU) is 3.4lnq % whereas the corresponding fraction for Sun-like 



stars is 7.6 ± 1.4% ( ICumming et al. 



20081 ). However, the sample size of planet-hosting M 



dwarfs is small and therefore the uncertainties are large. Further exploration of the giant 
planet population around M dwarfs may therefore offer some insight into which physical 
processes dominate formation of these worlds. 

Although the investigation of planets around M stars provides clues to their formation, 
the key to learning about how planets form is to observe them in the process. Recent 
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resu 


ts 


2005 


: ] 



t s on core accretion predict planet formation in mi 



Hubickyj et al 



2005 



Dodson-Robinson et al 



lions of years (e.g., lAlibert et al. 



2008). Gravitational instability models, 
however, predict ye ry fast f ormation times for massive planets in long period orbits: 10 3 



rtMayer et al. 



20041 ) to 10 5 flBodenheimer 



also occur on short timescales (~10 5 years, 



2006 ) years. Subsequen t inward migration can 



Papaloizou et al 



20071 1 that can bring these long 



period planets into the sensitivity range of RV surveys. Clearly, any data which can help to 
identify this timescale will help refine, limit, and possibly eliminate some planet formation 
models. Ongoing campaigns to catalog the planets around main sequence stars can not 
directly provide this information. By observing late-type, pre-main sequence (PMS) T 
Tauri stars, one acquires a snapshot of the early stages of planet formation around solar-like 
stars. 

Low-mass young stars present challenging targets for traditional RV surveys. Late 
spectral types, large distances (>100 pc), and extinction from natal dust clouds make 
these targets fain t at optical wavelengths. They also have strong magnetic fields (e.g., 



Johns- Krull 



20071 ) that generate large, cool star spots. These spots impact RV surveys of 



young stars by introducing significant Utter which can mimic the RV modulation imposed 



by a planet (ISaar fc Donahue 



1997 



Desort et al. 



20071 ). Recent attempts at detecting 
substellar compan ions in young stella r populations (10— 100 Myr) have generally been 



unsuccessful (e.g., 



Paulson et al. 



2004 : 



Paulson fe Yelda 



20061 ). likely because of the small 



sample size and intrinsic RV variability of the targets. The youngest RV planet detected 
to date is a ~6 M.tup planet on an 850 day orbit around the 100 Myr old star HD 70573 



(ISetiawan et al. 



20071 ) 



Spectral line bisector analysis has historically been used to disti nguish spot-induced 



RV variations from companion- induced ones (e.g.. lHatzes et al.lll997l ). Companion- induced 



reflex motion does not affect the shape of an absorption line's bisector whereas a spot 
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distorts the line symmetry. The line bisector span, the difference in bisector values at two 
different heights of the line profile, is a proxy for the average slope of the line bisector. 
A correlation between bisector span and RV variations suggests the RV fluctuation is 
spot-induced; otherwise, it may be caused by a companion. However, there may also be 
no correlation if the projected rotation velocity (v si n i) of the star is comparable to or 



less than the velocity resolution of the spectrograph (IDesort et al.l 



2007 



Prato et al. 



20081 ). 



Ther efore, bisector analys i s is only a first st ep in identifying potential young planet hosts 



e.g., 



Huelamo et al 



2008 



Prato et al 



20081 ) 



A potentially more reliable method for distinguishing between spots and planets 
leverages the wavelength dependence of the spot-induced RV modulation amplitude. The 
reflex motion induced by a planet affects all wavelengths equally. However, the contrast 
between a photosphere and a cooler star spot decreases at longer wavelengths because 
of the flux-temperature scaling in the Rayleigh- Jeans limit of blackbody radiation (e.g., 



Vrba et al. 



1986 



Carpenter et al 



20011 ). As a result, the amplitude of any spot-induced RV 
variability will be smaller at longer wavelengths. By observing in the visible and NIR it 
may be possible to distinguish between stellar activity and the presence of a companion by 



Reiners et al. 



(120101 ) point out that 



comparing the RV amplitudes at the two wavelengths, 
the magnitude of this decrease for spotted stars is dependent on the temperature difference 
between the photosphere and star spots and may not be significant if this difference is large 
(i.e. ~1000 K). Observations of T Tauri stars, however, do show reductio ns in RV jitter 



20081 ). lending support 



from optical to K band wavelengths by factors of 3 — 5 (jPrato et al.l 
to the plausibility of this approach. Furthermore, the late spectral types of T Tauri stars 
produce SEDs with peak emission around 1 /im and extinction of these partly embedded 



sources is lessened at longer wavelengths, thus allowing 



hence improved RV precision, in the NIR. 



Reiners et al. 



or in creased signal-to-noise, and 



( l2010l ) caution that the advantages 



of NIR RV measurements may not be apparent for stars earlier than ~M4— 5. For early 
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M stars, they assert that the decreased S/N in the optical bands is outweighed by the 
number of sharp spectral features available at the shorter wavelengths. Stars later than 
M5, however, do exhibit increased precision in the NIR both as a result of increased S/N in 
the J and H bands as well as the appearance of FeH absorption features. Finally, T Tauri 
stars allow for some relaxation in the desired RV precision. While optical surveys of main 
sequence stars strive for 1 m s " 1 precision or better, T Tauri stars exh ibit an intrinsic NIR 



RV variability of > 100 m s 1 (IPrato et al 



20081 : 



Mahmud et al. 



20111 ) . Therefore, NIR RV 



precision of several tens of meters per second is more than adequate for these targets. 
Thus, along with studies of planets around M and L dwarfs, RV surveys for young 



planets also drive the development of improved IN 



groups are active in th is area (e.g. 



2008 



Martin et al. 



IR RV precision measu r ements. Severa l 



2006 



Blake et al. 



2007 



Huelamo et al. 



Blake et al. 



the optical ( IBean et al 



2010 



with some reports of NIR velocity precision comparable to that in 



2010 



Figueira et al 



20101 ) . While these results are encouraging, 



these efforts are all focused on large (8 — 10 meter) aperture telescopes. There is therefore 
a need to develop precision NIR RV techniques for smaller aperture telescopes where more 
observing time is available to the community at cadences better tailored to RV surveys. 

In 2004, we began the M cDonald Observatory Young Planet Survey to monitor PMS 



stars in the nearby (~140 pc, 



Kenyon et al 



1994 ) Taurus-Auriga low-mass star forming 



region for evidence of substellar companions. Our sample of 143 c 



Tauri stars consists of V < 14 stars, most with v sinz < 20 km s 1 ( IHerbig fc Bell 



assical and weak-lined T 



19881 ) . For 



the first four years, a 



Coude Spectrometer (ITull et al. 



1 visible light o 



nervations were conducted using the Robert G. Tull 



1995|) on the 2.7-meter Harlan J. Smith Telescope. In 2008, 



we began to include ad ditional observations from the Sandiford Cassegrain Spectrograph 



( McCarthy et al 



19931 ) on the 2.1-meter Otto Struve Telescope. Observations at the Kitt 



Peak National Observatory (KPNO) 4-meter Mayall Telescope using the cassegrain echelle 
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spectrograph were also initiated in late 2008. Promising targets from these observations, 
based on measured RV variations and the results of bisector analysis, are selected for follow 
up K band observations at the 3-meter NASA Infrared Telescope Fa cility (1RTF) using th e 



high- resolution Ca ssegrain-mounted echelle spectrograph, CSHELL (ITokunaga et al. 



Greene et al. 



19901 : 



19931 ). Given the 1—3 Myr ages of our targets, a positive detection of a young 



exoplanet will provide a unique datapoint for giant planet formation theory. 

In this paper, we discuss our methodology for measuring precise NIR RVs with 
CSHELL. In §[2] we discuss our observing strategy and data reduction algorithms. In §0 
we present our methodology for using telluric absorption features to measure the RVs of 
our targets. We present results from observations of an RV standard star and a known 
exoplanet in and we discuss the limitations and possibilities of our technique in $5j 



2. Observation strategy and data reduction 



Observation s were taken at the 3 meter NASA Infra red Telescope Facility (IRTF) 



using CSHELL (ITokunaga et al. 



1990 



Greene et al. 



1993|). CSHELL is a long-slit echelle 



spectrograph (1.08 /im — 5.5 fim) that uses a Circular Variable Filter (CVF) to isolate a 
single order onto a 256x256 InSb detector array. We used the CVF to isolate a 50 A segment 
of spectrum centered at 2.298 /im. This region contains numerous deep photospheric 
absorption lines from the CO v = 2—0 bandhead as well a rich set of predominately CH 4 
telluric absorption features which we use as a wavelength reference. The 0".5 slit yielded a 
typical FWHM of 2.6 pixels (~0.5 A, measured from arc lamp spectra) corresponding to a 
spectral resolving power of R ~ 46,000. 



We observed RV standards, a known exoplanet candidate, and several T Tauri planet 
host candidates in February 2008 (8 nights), November 2008 (6 nights), February 2009 (2 
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nights), November 2009 (5 nights), and February 2010 (8 nights). At the beginning of 
each night, we imaged twenty flat fields, each with a 20 second integration time, using a 
continuum lamp to illuminate the entire slit. We also imaged the same number of 20-second 
dark frames. Additionally, we imaged six Ar-Kr-Xe emission lines, changing the CVF while 
maintaining the grating position, to determine the wavelength reference. All of our target 
data were obtained using 10" nodded pairs to enable subtraction of sky emission, dark 
current, and detector bias. Integration times for each nod were typically 600 seconds; for 
fainter targets we took multiple contiguous nod pairs. The signal-to- noise ratio (S/N) for 
all targets varied significantly depending on cloud cover, seeing, guiding errors, etc., with 
typical values for all targets of ~50/pixel per nod position. 



Johns-Krull et al. 



(1999) 



Our data reduction strategy closely follows that described in 
and was implemented entirely in IDL. We first produced a nightly master dark frame by 
averaging the twenty individual dark exposures and found a median dark current of 0.24 
e~/s. We produced a nightly normalized flat field by averaging the flat field exposures, 
subtracting the master dark frame, and then dividing by the mean of the dark-subtracted 
master flat. For our target data, we first subtracted the nod pairs to create a difference 
image. The difference image was subsequently divided by the normalized flat field. We 
estimated the read noise from the standard deviation of a Gaussian fit to a histogram of 
the pixel values in the difference image (~30 e~). We then identified the location of the 
curved spectral traces in the difference image by fitting a second-order polynomial to the 
location of maximum and minimum flux along the dispersion direction. 

To optimally extract the spectrum, we divided each nod pair into four equally spaced 
bins of 64 columns along the dispersion direction. Within each of these 64 column bins 
we constructed a lOx oversampled "slit function" (i.e., the distribution of flux in the 
cross-dispersion direction). First, a rough estimate of the spectrum was created by summing 
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the pixels in each column of the difference image for each nod position. The limits included 
in this sum are from the midpoint between the two nod positions on the detector to the 
edge of the area that is well illuminated by the flat lamp. This is typically 60-70 pixels in 
each column for each nod position. Next, each pixel in the bin was sorted by its distance 
from the order center for the column the given pixel falls in. The flux in each pixel was 
divided by the rough estimate of the spectrum for its appropriate column to normalize all 
the pixels going into the slit function estimate. The flux in these offset ordered pixels was 
then median filtered with a 7 point moving box. A flux estimate for each oversampled pixel 
was then determined by taking the median of all the pixels that fell in a given subpixel. 
This then formed our oversampled master slit function. The multiple median filters 
generally remove the effects of cosmic rays and uncorrected bad pixels on the determination 
of the slit function. We then fit this master slit function to a three Gaussian model: a 
central Gaussian flanked by two satellite Gaussians. The amplitude, center, and width of 
each Gaussian were fit as free parameters using th e IDL implementation of the AMOEBA 



non-linear least-squares (NLLS) fitting algorithm (INelder fc Meadl 119651 ) . The resulting 
model was then normalized to unit area. This algorithm produces four model slit functions, 
one for each bin. However, the actual slit function is a smoothly varying function of column 
number. Therefore, to smooth out the transitions from bin to bin, we create 256 column 
slit functions by linearly interpolating between the four bin slit functions. 

To determine the total flux in each column of the spectrum, we calculated the scale 
factor that best m a tches our model slit function to the column data, per the recipe 



described in iHornd ( 119861 ). In order to mask out spurious flux levels from cosmic rays, we 
implemented an iterative sigma-clipping algorithm. This algorithm starts with an estimate 
of the total noise from the measured read noise in the differenced image plus the Poisson 
noise from the target. We then subtracted our model fit from the data in each detector 
column and masked those pixels whose residual was 3-a greater than the estimated noise. 
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We iterated through this process 1—2 times until the scale factor converged, thus providing 
an optimal value of the spectrum in that column that is largely immune to hot pixels, 
cosmic rays, etc. This algorithm also provides an estimate of the flux uncertainty at each 
location along the spectrum. 



3. Radial velocity determination 



The telluric absorption features in the K-band provide an absolute wavelength and 
instrumental profile reference, similar in con cept to the iodine gas cell technique used in 



high-precision optical RV exoplanet surveys (IButler et al. 



19961 1 . Using the atmosphere as 



a "gas cell" lets us superimpose onto our spectra a relatively stable wavelength reference 
which follows the same optical path as the light from the science target. This helps alleviate 
uncertainties introduced by variable slit illumination, changing optical path lengths, etc. 
We determined the radial velo cities of our target s using a spectral m odeli ng technique 



simila r to the one presented in 



Blake et al. 



(120071 ). 



Prato et al. 



20081 1 . and 



Figueira et al. 



( 120101 ). We modeled the stellar spectrum and the telluric features using two high-resolution 



template spectra. For the stel lar spectrum, we employed NextGen stellar atmosphere 



models (IHauschildt et a 



We used SYNTHMAG 



1999) t ailore d to the T e ft, logg, and metallicity of our targets. 



Piskunov 



with atomic ( IKupka et al.l 120001 ) and CO ( IGoorvitch 



19991 ) to generate spectra from the NextGen models along 



19941 ) line lists. These spectra are 



sampled at a resolution ~14 times greater than our obs ervations. We modeled the telluric 



features using the NOAO telluric absorption spectrum ([Livingston fc Wallace 



199J) which 



has a resolution ~4 times higher than that of our observations. 

In order to match the model to our observations, we fit a velocity shift and a power law 
scaling factor for each template separately. The power law scalin g accounts for difference s 



in line strength between the observed and template spectra (e.g., 



McCullough et al. 



20061 ) 
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The stellar spectrum was then interpolated onto the wavelength scale for the telluric 
spectrum. To match stars with measurable rotation, we fit for the value of vsini that 
optimally broadens the stellar template spectrum using a disk-integration routine with 
a limb-darkening coefficient of 0.2 (based on fitting the model photosphere to a linear 
limb-darkening law). For the main sequence stars observed for this paper, however, 
we fixed vsini to zero since it is not meaningfully detected above the instrumental 
broadening. Because the normalized telluric spectrum models the relative transmission of 
the atmosphere, we multiplied the telluric spectrum and broadened stellar spectrum together 
to create a composite spectrum and convolved the result with a Gaussian instrumental 



profile (IP) to model t 
models of the IP (e.g., 



re instrument response. We experim ented with multi- Gaussian 



Valenti et al. 



1995 



Bean et al. 



20071 ) but found no improvement in 



the RV precision. We then fit the continuum level with a second-order polynomial. Finally, 
we binned the 4x oversampled composite model down to the resolution of our CSHELL 
data. The resultant model has nine free parameters (two velocity shifts, two scaling factors, 
vsini, IP FWHM, and three continuum coefficients). To simultaneously determine the best 
fit to these parameters, we again used the IDL implementation of the AMOEBA NLLS 
algorithm. 

Because CSHELL is cassegrain-mounted, instrument flexure results in a wavelength 
dispersion on the detector that changes with each observation. We determined an initial 
dispersion solution using the Ar-Kr-Xe lamp lines imaged at the start of each night (§2J). 
The six individual arc lamp spectral images were stacked to create a composite spectrum. 
The combined arc spectrum was then extracted using the order tracing from that night's 
highest S/N stellar spectrum. We found an initial dispersion solution by fitting a third 
order polynomial to the locations of these emission lines. This initial polynomial serves as 
a starting point to refine the dispersion solution. For each target observation, we take the 
observed spectrum and the best-fit model and divide both into eight, 6.25 A wide, bins. 



13 



This breaks up the observed and model spectra in to eight pieces. We then take each piece 
of observed spectrum and cross-correlate it with the corresponding piece of model spectrum 
to measure the pixel shift between the two. The pixel shift as a function of spectrum bin 
provides a crude estimate of how much the observed dispersion solution differs from the 
initial solution across the detector. The eight pixel shift values are fit to a second-order 
polynomial. By interpolating this polynomial across all detector columns, we estimate by 
how much the wavelength has shifted per pixel relative to the initial dispersion. A new 
wavelength solution is determined by adding this shift to the initial solution. We then 
produce a new model fit using the refined wavelength solution. We iterate this algorithm 
~5 times until the wavelength solution converges. 

Figure [1] illustrates the outcome of this modeling for one observation of our RV 
standard, GJ 281. Radial velocities for each observation were calculated by taking the 
difference in the velocity shifts between the best-fit stellar and telluric spectra. These 
velocities were in turn corrected for barycentric motion determined for the mid-time of each 



exposure Q 

For each observation, we estimated uncertainties introduced by the photon-limited 
errors of our observations. These uncertainties in turn have two main components: the 
first is from the information content in the stellar spectrum itself and the second is in the 
information content in the tellu ric spectrum which is used as a wavelength calibrator. As 



described in 



Butler et al 



(jl996[ ) , the intrinsic RV error of a spectrum is the quadrature sum 
of these two error contributions. This in turn is a function of the slope of the intensity with 
respect to wavelength at each pixel and the S/N. For a given S/N, a spectrum with more 
numerous and sharper absorption features will have more information content than one with 



2 Barycentric corrections were calculated using helcorr.pro, an IDL routine based on 
the IRAF task noao . astutils . rvcorrect. 
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fewer, broader features. We used the best-fit model photosphere to calculate the intensity 
derivative and the measured flux uncertainty (estimated from the measured read noise plus 
the target Poisson noise, §2]) to determine the S/N for each pixel. The combined errors from 
each pixel determine the photon-limited uncertainty on the final RV. We derived the errors 
introduced by the telluric spectrum in exactly the same manner, replacing the photosphere 
model with the telluric one. The final uncertainty on the velocity was calculated by adding 
the stellar and telluric errors in quadrature. 

A final, nightly RV was determined by calculating the average of the individual nod 
RVs, weighted by the median S/N across each spectrum. The uncertainties were computed 
by taking the weighted standard deviation of the nod RVs and dividing by the square root 
of the number of nods. 



4. Results 



4.1. GJ 281 



To determine the precision of our approach, we obtained 28 observations of GJ 
28 1 over the past two yea rs. This is a late-type star (MO, K mag = 5.9, T^ff = 3776 



K. 



Casagrande et al. 



20081 1 known to have a stable RV (r.m.s ~ 12 m s , 



Endl et al. 



20031 ). We determined the nightly RVs and uncertainties using the approach described in 
Section [3] (see Figure [2] and Table [I]). The median uncertainty from photon statistics is 
40 m s _1 per nod with the telluric uncertainty ~50% greater than the stellar spectrum 
uncertainty. The standard deviation of the mean within a sequence of contiguous nod pairs 
{&Rv/^fN) is typically 37 m s" 1 . After averaging, the standard deviation of the nightly 
RVs over the entire 24 month observing window is ~58 m s -1 . Since this is greater than 
the variability seen in a given night, we need to estimate what additional systematic error 
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we are introducing between nights. Assuming that the night-to-night systematics add in 
quadrature with the variability over multiple contiguous nods to give us our overall 58 
m s" 1 standard deviation, {<J% OT = o 2 nods + cr 2 ight i y ), we estimate y/58 2 — 37 2 = 45 m s _1 as 
our night-to-night systematic limit. 

To better understand the systematic errors, we looked for differences in the measured 
velocities as a function of nod position. Subtracting the A beam velocities from the B 
beam velocities, we found a mean difference between the beams of 66 m s -1 . The standard 
deviation of these differences was 131 m s _1 , while the standard deviation of the mean was 
22 m s^ 1 . The standard deviation of the mean of just the A beam velocities was 15 m s _1 ; 
for the B beam velocities, this was 17 m s _1 . This suggests that there is a statistically 
significant difference between the velocities at the two nod positions which gets averaged 
out at some level when combining all observations in a single night. A simple subtraction 
from the A velocities to remedy this offset is not sufficient, however. In a few nod pairs, 
the B velocities are higher than the A velocities and differencing makes this offset worse. 
Because of this, the variability seen in the binned velocities is nearly identical whether or 
not the mean offset is removed. 

It is not clear why there is a velocity offset between the beams. While we expect to see 
different distortions in the two nod positions, the use of telluric features as a common-path 
calibrator should account for this. One concern with our observing setup is the asymmetric 
distribution of telluric lines across the spectrum. The long wavelength half of our spectra is 
much richer in telluric features than the short- wavelength half. It is therefore possible that 
we are creating false higher order terms as we track the change in the dispersion across the 
slit. We tested this by recalculating our RV determination for all the GJ 281 observations 
using only the left half of the spectrum and then using only the right half of the spectrum. 
Fitting only the left half results in a mean difference between the A and B beams of 33 
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m s _1 whereas restricting ourselves to the telluric-rich right half results in a mean difference 
between the beams of 10~ 6 m s" 1 . While this would seem to support the hypothesis that 
an asymmetric telluric line distribution may be introducing systematic errors, the standard 
deviation of the mean in each beam for both cases is large enough to make the statistical 
significance of these results unimportant. Interpretation of the much lower RV scatter 
in the right half of the spectrum is therefore unclear. It is likely that we simply do not 
have enough signal to noise to effectively separate out this error from our other sources of 
uncertainty. Future observations of very bright targets may help identify the sources of our 
systematic error. We note that the existence of a systematic shift between the two beams 
combined with our use of the standard deviation of the measured RVs in both nod positions 
will most likely over estimate the uncertainty in our velocity measurements. 



4.2. Gl 86 



In addition to RV standards and T Tauri planet host candidates, we also obtained six 
observations of the kn own exoplanet host Gl 86 (K mag = 4.1 , T e fj = 5350 K, log g = 4.6, 



Flvnn & Morel! 



19971 ) in November 2008. 



Queloz et al. 



( 120001 ) prese nted evidence for a ho t 



Figueira et al. 



Jupiter around Gl 86 with msini = 4.0 Mjup and P = 15.78 days, 
used Gl 86 as a test case for their precision NIR RV work with the CRIRES spectrograph 
on the VLT and were able to show that their observations were consistent with the 



Queloz et al. 



orbital solution, allowing for errors in time of periastron (To) determination 



and drift in the center-of-mass veloc ity (V r ) as a result of the presence of a long period 



companion ( lEggenberger et al. 



20031 ). We present our nightly RVs for Gl 86 in Table M 



and Figure El The median uncertainty from photon statistics is 41 m s _1 with the telluric 
uncertainty ~43% greater than the stellar spectrum uncertainty. The standard deviation of 
the mean within a sequence of contiguous nod pairs is typically 58 m s -1 . Figure |3] plots our 
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binned RVs with those from iFigueira et all and an orbital fit to both datasets. Our error 
bars are determined by adding the results from our binning algorithm (§3J) in quadrature 
with the 45 m s _1 systematic uncertainty determined from our RV standard ( §4.ip . When 
fitting an orbit model to the data, we allowed T n a nd V r to be fit as free parameters, while 



fixing all other values to those from 



Queloz et al 



((20001) • We show good agreement with 



previously published values exhibiting a ~46 m s 1 standard deviation in the residuals. 

All of our RVs differ from the orbital solution by no more than 0.5 sigma; the 
mean difference is ~0.24 sigma. Assuming normally distributed errors, we would expect 
1 — 2 of these residuals to be greater than one sigma. To explore this further, we 
assumed that the underlying uncertainties are gaussian and drawn from a distribution 



with cr = V58 2 + 45 2 = 73 m s _1 . The variance in the residuals is dominated by the first 
RV measurement; excluding that point we find a standard deviation in the residuals of 
16 m s -1 . We used a Monte Carlo analysis to estimate that the probability of having at 
least five out of six points within 16 m s _1 of the mean, given a = 73 m s _1 , is ~0.08%. 
We further explored the false alarm probability (FAP) that the phase coherence from our 
observations is the result of chance. We performed a Monte Carlo analysis by holding the 
times fixed and sampling the CSHELL velocities (and corresponding uncertainties) with 
replacement over 10 4 iterations. For each iteration, we fit an orbit model as described above 
and recorded how often the reduced \ 2 of the fit was less than that of the fit to the observed 
velocities. We find that it is highly unlikely that the phase coherence is a result of chance, 
with a FAP of 0.0003. We are therefore very confident that our RVs are consistent with the 
known planetary-mass object orbiting Gl 86. The inconsistency between our estimated error 
bars and the residual scatter is suspicious, however. The scatter in the residuals is more 
comparable to the estimated photon uncertainty. Our 45 m s -1 systematic uncertainty is 
measured over a period of two years whereas the RVs for Gl 86 were measured over a single 
run. If the systematics were "quiet" on this run, then 45 m s _1 may not be the appropriate 
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value to use. We see some evidence for the appearance of "quiet" times in our data on GJ 
281. During the November 2008 run the standard deviation of our nightly (6 nights) GJ 
281 data is 39 m s _1 , well below the 58 m s~ x of the entire data set. This is the same run on 
which the Gl 86 data were obtained. We therefore believe that our uncertainties are most 
likely overestimates. 



5. Discussion 



NIR spe ctroscopy has only r ecently begun to play a role in the search for substellar 



companions. 



Martin et al. 



(120061 ) combined optical and H band observations to look for 



giant planets around the young brown dwarf LP 944-20; they concluded that the observed 
optical RV modu 



Blake et al 



ations were driven by inhomogeneous surface features (i.e. clouds). 
(120071 ) used high resolution K-band spectroscopy to investigate the presence of 
giant planets around a population of L dwarfs, achieving a precision of 300 m s" 1 . They 
find no evide nce for companions w ith Msinz > 2Mj and P < 3 days in their sample of 



nine targets. 



Setiawan et al. 



(120081 ) reported the detection of a ~10 Mj planet on a 3.5 



da y orbit around the 10 Myr old star TW Hydra. However, H band observations presented 



by 



Huelamo et al. 



(120081 ) reveal a strong wavelength dependence in the velocity amplitude, 



thus casting doubt on the presence of a compan i on an d suggesting that spots are the 



cause of TW Hydra's RV variations. iPrato et al 



(120081 ) observed two potential planet-host 



candidates in the NIR, selected on the basis of their optical variability. To within their 
measurement precision, they detect no K band RV modulation, implying that spots are 
responsible for the apparent optical RV variability. 

The aforementioned interest in targeting M dwarfs as hosts of habitable planets ( JQ) 
has spurred yast i mprovements in NIR precision, bringing it closer to that of optical surveys. 

(120101 ) report achieving a long term precision in the H band of ~5 m s _1 on late 



Bean et al 
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M dw arfs using CRIRES at the VLT with the aid of an ammonia gas cell. iFigueira et al 



( I2010l ) report a comparable precision also on CRIRES but using t elluric absorption features 



Blake et al 



(120 id ) report 



as a wavelength reference. Based on six years of NIRSPEC data, 
a precision of 50 m s" 1 for their sample of K dwarfs and 200 m s _1 for L dwarfs also through 
the use of telluric lines. 



Prato et al. 



fcoosj ) 



In an earlier paper from our survey of young, low-mass stars, 
presented evidence that the RV variability of two T Tauri stars, DN Tau and V836 Tau, is 
the result of star spots. These two targets exhibited an optical RV standard deviation of 
438 m s" 1 and 742 m s _1 , respectively. K band RVs, however, showed a standard deviation 
of 144 m s _1 for DN Tau and 149 m s -1 for V836 Tau. These values are consistent with 
the measurement uncertainties of ~165 m s _1 and ~290 m s _1 , respectively. These stars 
therefore exhibit no evidence of RV jitter in the K band. The undetected NIR RV variability 
demonstrates two important points. First, optical RV surveys can not rely on bisector 
analysis alone to ascertain the mechanism for RV variability. Both DN Tau and V836 Tau 
showed significant periodic optical variability and neither exhibited a correlation between 
RV and bisector span. Multiwavelength observations are therefore an essential element for 
testing the planet hypothesis, especially in young stars. Second, they support the hypothesis 
that RV jitter from stellar activity should be significantly lower at longer wavelengths. We 
have demonstrated ~58 m s _1 precision with CSHELL over two years and the ability to 
detect a known 4 Mjup planet. If the ~150 m s -1 variability observed in DN Tau and 
V836 Tau is typical of T Tauris, our 58 ms" 1 precision is well below the intrinsic noise floor 
of our targets suggesting that our methodology is well-suited to detecting "hot Jupiters" 
around young stars. However, we caution against extrapo l ating from two young stars to 



the NIR behavior of other T Tauri stars. 



Mahmud et al. 



(120 111 ) present another target 



from our Taurus survey that exhibits a ~430 m s 1 modulation in the K band, which we 
believe to be caused by spots, suggesting that the NIR RV behavior of T Tauris may vary 
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significantly from star to star. 
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Fig. 1. — Illustration of NIR spectrum modeling for our highest S/N observation of GJ 
281 (JD 2455252.906): (a) NextGen photosphere model, (b) telluric template, (c) combined 
model (solid line) with CSHELL observation (dots), (d) residuals. All plots are normalized 
to the final model and have constants added for visual clarity. 
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Feb 2008 Nov 2008 Feb 2009 Nov 2009 Feb 2010 




Days since Feb 1 3, 2008 

Fig. 2.— Binned K band radial velocities of GJ 281 from Feb 2008 - Feb 2010. The error 
bars are the standard deviation of the mean of all nods obtained on that night (<J no d s / VN)', 
they exhibit a median value of ~37 m s _1 . The standard deviation of the binned RVs over 
the entire observing period is 58 m s" 1 . 
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ig. 3. — Known exop lanet Gl 86 b. Open blue squares are NIR RVs measured using CRIRES 



Figueira et al 



2010|). Filled red circles are our binned K band RVs from November 2008. 



Black line is orbit model fit to the combined datasets. Weighted standard deviation of our 
residuals is 46 m s _1 . The bar in the lower left indicates the median uncertainty from photon 
statistics. 
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Table 1. GJ 281 Radial Velocities 



JD - 2450000 Heliocentric RV (km s" 1 ) a nods /VN (m s" 1 ) 



4510.914 


20 


,096 


70.8 


4511.910 


19. 


.985 


14.2 


4512.914 


19 


.981 


67.7 


4513.895 


20 


.089 


37.4 


4514.910 


20 


.092 


60.6 


4515.902 


20 


.054 


5.5 


4516.902 


20 


.129 


51.6 


4781.137 


20 


.108 


12.3 


4783.125 


20 


.134 


30.7 


4784.109 


20 


.161 


26.0 


4785.082 


20 


.114 


38.4 


4786.109 


20 


,048 


71.0 


4787.133 


20 


.138 


59.0 


4884.906 


20 


.065 


28.0 


4885.867 


20 


.126 


29.3 


5151.137 


20 


.185 


41.8 


5153.160 


20 


.137 


50.9 


5156.148 


20 


.094 


6.4 


5158.152 


20 


.054 


30.7 


5160.145 


20 


.161 


19.0 


5235.934 


19 


.988 


42.7 
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Table 1 — Continued 



JD - 2450000 Heliocentric RV (km s^ 1 ) a nods /VN (m s" 1 ) 



5236.922 19.996 32.8 

5237.922 20.171 26.3 

5238.926 20.135 10.8 

5239.918 20.053 61.0 

5240.926 20.150 63.8 

5241.906 20.107 12.8 

5242.902 20.066 86.6 
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Table 2. Gl 86 Radial Velocities 



JD - 2450000 Heliocentric RV (km s" 1 ) a nods /VN (m s" 1 ) 



4781.902 56.303 199.5 

4782.891 56.087 73.0 

4783.879 56.008 52.8 

4784.883 55.974 52.3 

4785.891 56.054 65.4 

4786.891 56.118 81.2 



